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Chapter 1 
Introduction 


Thin sheet flows have been a matter of research interest for many years. Sir 
Geoffrey Taylor has published a series of papers on the dynamics and stability of liquid 
sheets, [1]. Dombrowski and Fraser have written a paper on the disintegration of 
expanding cylindrical sheets, such as would be found in a fuel injector, [2], Related work 
has been done by Crapper on capillary waves on fluids of infinite depth, [3]. As long ago 
as 1879, Rayleigh published a paper on the capillary break-up of a cylindrical liquid jet in 
air, [4]. Interest in thin sheet flows has recently been renewed due to their potential 
application in space radiators. 

1.1 - Historical Perspective 

A classic illustration of the methods of hydrodynamic stability is Rayleigh's theory 
of the break-up of a liquid jet in air, such as the formation of drops by a thin jet of water 
from a faucet. It had previously been shown that capillarity would lead to instability of a 
round jet, because an axisymmetric deformation could decrease the surface area of the jet . 
Rayleigh analyzed in detail the instability of a uniform basic flow of incompressible, 
inviscid liquid within a cylinder; the liquid having a free surface governed by its surface 
tension. The jet was found to be stable to all non-axisymmetric modes, but is unstable to 
axisymmetric modes whose wavelength is greater than the circumference of the jet 
Rayleigh's theory agrees well with experiments. Rayleigh also found that there is a 
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recognized word, varicosity, to describe the axisyrametric instability. This is sometimes 
called sausage instability now, [4J. 

In the second paper of Geoffrey Taylor’s series he discussed the waves on thin 
sheets of fluid. He showed that capillary waves are of two kinds, symmetrical waves in 
which the displacements of opposite sides are in opposite directions, and antisymmetrical 
waves in which the displacements are in the sajne direction, as shown by Figure 1-1. Any 
disturbance can be regarded as composed of these two types of waves. The 
antisymmetrical waves were shown to be non-dispersive. The antisymmetrical waves in a 
sheet of uniform thickness caused by a point disturbance were shown to appear as two 
narrow line-like waves. The symmetrical waves were very different, they were highly 
dispersive and were propagated much more slowly than the antisymmetrical waves. 
Experimentally, a point disturbance produces both kinds of waves simultaneously. The 



Symmetric Waves 



Asymmetric Waves 

Figure 1-1 Symmetric and Asymmetric Waves 
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symmetrical waves in a sheet of uniform thickness caused by a point disturbance were 
shown to appear as parabolas, [1], 

In 1960, D. R. Brown carried out an experimental investigation of thin liquid sheet 
flow in connection with a method of lacquer application known as curtain coating. In 
curtain coating, the coating material passes through a thin slit and a liquid sheet is 
produced. The width of the sheet is maintained by guide wires along the edges of the 
sheet. The minimum liquid flow rate required to maintain a stable sheet is discussed, and 
the effects of the impingement of the sheet on a rapidly moving surface are described, [5]. 

The stability of a viscous liquid curtain falling down steadily under the influence of 
gravity was investigated by S. P. Lin in 1980. A linear stability analysis was carried out 
and it was found that only the spatially growing disturbances whose group velocity points 
toward the upstream direction were unstable. A critical Weber number was found. Any 
sheet formed with a Weber number less than the critical Weber number was found to be 
unstable, [6], This was in good agreement with the experimental results found by Brown. 

1.2 - Space Application 

One way of significantly reducing the mass of a space radiator is to eliminate the 
containing walls for the working fluid. Thus the working fluid is exposed to the vacuum 
condition of space. Such an external flow radiator will have a lower mass than a heat pipe 
type radiator. Ease of deployment and near immunity to meteoroid damage are two other 
advantages of external flow radiators. 

Currently there are three different designs for external flow radiators. They are the 
liquid droplet radiator (LDR), the liquid belt radiator (LBR), and the liquid sheet radiator 
(LSR). The liquid sheet radiator uses a thin liquid sheet as the radiating surface. Tne sheet 
may be as thin as 0.0015 in. 


.S' 


The idea behind the design of the liquid sheet radiator is to have the greatest surface 
area for the least volume or mass. With a very thin liquid sheet the mass is very small, 
while the surface area, both sides of the sheet, is quite large, allowing for much more 
specific surface area than the LDR or LBR. Also, being a single sheet, the radiation of heat 
from one droplet to another droplet is eliminated. 

One of the advantages of the LSR is ease of design. The machining required for the 
fabrication of the narrow slits that are used in producing the sheet does not require as much 
precision as that for the fabrication of the many small holes for the liquid droplet radiator. 
Also, the pump power will be less since the viscous losses for a single slit are less than that 
for many small holes. Collection of the liquid is also greatly simplified since the sheet 
coalesces to a single point, as will be shown later. The simpler design for the sheet 
generator and collector should result in a lower mass for these components, [7J. 

Due to its many advantages the LSR has been the subject of great interest in recent 
years. This has resulted in the need to fully understand the fluid dynamics and stability of a 
thin liquid sheet. 

1.3 - Introduction to Liquid Sheet Flow 

A liquid sheet is formed by forcing the working fluid through a very thin sliL Since 
the sheet is so thin, the surface tension has a very pronounced effect. The sheet edges are 
pulled towards the center of the sheet, (see Figure 1-2). This causes the sheet to take on a 
triangular form. Due to conservation of mass the fluid must collect in the so-called edge 
cylinders that border the sides of the sheet. 

In order to more fully understand the fluid dynamics of liquid sheet flow, a number 
of analytical and experimental studies were carried out The first analytical study, Chapter 
2 Section 1, was to determine an analytical prediction for the sheet length over slit width 
ratio for the liquid sheet. The second. Chapter 2 Section 2, was to determine the theoretical 
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Figure 1-2 Photograph of a Thin Liquid Sheet 

T cross-sectional shape of the edge cylinders. 

In addition to the fluid dynamics of liquid sheet flow it is also of vital importance to 


i 


y- understand the Suability of the flow, i nis is the focus of Chapter 3. Tn Chapter 3, a linear 
Jjjk'\ stability analysis is presented on a nor. -planer liquid jet. This analysis shows the possible 
instability of the edge cylinders. 





Chapter 4 has two main parts. In the first part, the experimental apparatus is 
described. In the second part, the formation of holes in the sheet is discussed and the 
method of eliminating those holes is presented. 

Chapters 5 and 6 present the experimental procedures and the experimental results, 
respectively. Experiments were carried ounp measure the length over width ratios for the 
liquid sheet and comparisons were made to the theoretical predictions presented in Chapter 
2. An elaborate photographic technique was employed in order to determine the cross- 
sectional shape of the edge cylinders. This experiment revealed a very interesting 
phenomenon, the edge cylinder cross-sectional shape oscillates between an elliptical and a 
peanut-like shape, and the theoretical predictions are only valid in the limiting case. A third 
experimental study was performed on the effects of air drag on a liquid sheet. 



1.4 - Introduction to Linear Stability Analysis 

The essential problems of hydrodynamic stability were recognized and formulated 
in the nineteenth century, notably by Helmholtz, Kelvin, Rayle.igh and Reynolds. The 
method of normal modes for studying the oscillations and instability of a dynamic system 
of particles and rigid bodies was already highly developed. A known solution of Newton's 
or Lagrange’s equation of motion for the system was perturbed. The equations were 
linearized by neglecting products of the perturbations. It was further assumed that the 
perturbation of each quantity could be resolved into independent components or modes 
varying with time t, like e st for some constant s, which is in general complex. The values 
of s for the modes were calculated from the linearized equations. If the real part of s was 
found to be positive for any mode, the system was deemed unstable because a general 
initial small perturbation of the system would grow exponentially in time. 

Stokes, Kelvin, and Rayleigh adapted this method of normal modes to fluid 
dynamics. The essential mathematical difference between fluid and particle dynamics is 
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that the equations of motion are partial rather than ordinary differential equations. This 
difference leads to many mathematical difficulties in hydrodynamic stability, which have 
been overcome for only a few classes of flow with very simple configurations. 
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Chapter 2 

Theoretical Formulation 

The focus of this chapter is to analytically determine the sheet length over slit width 
ratio and the cross-sectional shape of the edge cylinders. A better understanding of the 
problem at hand can be realized by referring to Figure 2-1, where the front and cross- 
sectional views of the liquid sheet are shown schematically. The sheet is formed by forcing 
the working fluid through a very thin slit, of length Wq and thickness z q. The sheet is then 
pulled together by surface tension until the sheet coalesces at a length L. An orthogonal- 
cartesian coordinate system centered at the top of the sheet in the middle of the slit is 
adopted. The z-axis points in the vertical down direction, the x-axis points to the right 
along the slit, and the y-axis goes in the direction of the thickness of the sheen The 
velocity convention that u, v, and w denote the x, y, and z components of velocity, 
respectively, will be used throughout the chapter. Moreover, s(x,z) represents the function 
that will describe the cross-sectional shape of the edge cylinders. Using the above 
mentioned notation, the two objectives of this chapter are to determine L/W 0 and s(x,z). 

In the following section the effects of gravity and surface tension on the sheet will 
be considered. First, the forces acting on the edge cylinder will be considered, resulting in 
an expression for the velocity of the edge in the x-direction. Second, a macroscopic mass 
balance on the entire sheet will be undertaken to determine the sheet length over slit width 
ratio. Third, the sheet length over slit width ratio will be found by neglecting gravity and 
compared to the ratio that included the effects of gravity. This analysis will result in 
showing that the surface tension is the dominating force in the sheet and that the 
gravitational force may be neglected. 
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In the second section the only force acting on the sheet will be the surface tension 
force. Here, the equations of continuity and motion will be solved to obtain the cross- 
sectional shape of the edge cylinders. A number of assumptions will be made resulting in 
the motion equation just being a balance of inertial forces and surface tension forces. The 
objective is to develop a method for predicting the shape of the edge cylinder from knowing 
the area of the edge cylinder. This area, in turn, _ can be found from the results of the first 
section. 
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Front View 


Cross-Sectional View 


Figure 2-1 Schematic of Thin Sheet Flow 


2.1 - Derivation of the Sheet Length over Slit Width Ratio 


2.1.1 - Development of the Governing Equations 


Consider an arbitrary element of the edge cylinder with a thickness of dz, as shown 
in Figure 2-2. Let this element be considered a control volume. Assume it to be at pseudo- 
steady state, that is, it will be considered steady when conservation of momentum is 
applied and unsteady when conservation of mass is applied. Furthermore assume it to be 
unaffected by gravity. Assuming a zero slope at the point where the cylinder connects to 
the sheet, the surface tension force between the cylinder and the sheet is adz for both the 
front and the back of the sheet, therefore: 

dF x = 2 adz 


or 


dz 



Figure 2-2 Differential Element of the Edge Cylinder 



Let be equal to the velocity at which fluid enters the control volume at the right 
edge. Applying a macroscopic momentum balance on this control volume yields: 

dF x = pu£xdz 


on 


Therefore: 


p *-^- 20 


u .-„/2a 


( 2 - 1 ) 


If Ae is defined as the area in the xy plane of the control volume, then applying a 
macroscopic mass balance on the control volume yields. 


- 4 - (pAcdz) = pivtdz 

dt 


Therefore: 




dt V p (2-2) 

If the principle of conservadon of mass is applied to the entire sheet (see Figure 2- 
1) then it can be shown that: 


p wo \Vo to — pw(2Ac + Wx) 

Since the fid d is incompressible and assuming to - W above equation can be 

written as: 


wo Wo to = 2 wAc + w 0 t 0 W (2-3) 

Since the change with unw of the sheet width is twice the edge velocity it may be 
written that: ~— 

2U' = ‘J* sl 

<* W px (2 . 4) 

It is known that the accelerate* in ** is due to gravity on 
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Integrating, knowing that at t = 0; dz/dt = w 0 , gives: 


^ = w = gt + w 0 

Taking equation (2-3) and differentiating with respect to time yields: 
0 =2 ^s. w + 2A. + Town ^ 

Differentiating again yields: 


d. ■■ ■“•f+w-dt 


( 2 - 6 ) 


But, 


•Resulting in: 


O = 2^w + 4^^ + 2A c ^ + t 0 w 0 ^ 
dt 2 dt dt ^ dt 2 dt 2 


d—w _ d§ _ a 
dt 2 dt 


0 = 2^w + 4^dw +t £W 

dt 2 dt dt dt 2 


Differentiating equations (2-2) and (2-4) with respect to time gives: 


d 2 Ac 

/ o dx. 

dt 2 * 

V 2px dt 

£SL m . 

/ .22.sk 

dt 2 ^ 

/ pt 3 dt 


(2-7) 


( 2 - 8 ) 


(2-9) 


’^Substituting equations (2-2), (2-5), (2-6), (2-8), and (2-9) into equation (2-7) yields: 

% 



0 = 2 V& * <gt + wo) + V 2 ? 8 + t ° wo 'V^ r * 

^Solving this for dx/dt gives: 

•4T 2 g 


di = _ 
dt x (gt + w 0 ) + Tqwo 


(2-10) 
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From this expression it is apparent that if g = 0 the thickness is constant, as would be the 
case in space application. But, since g is greater than 0, dx/dt is less than 0 and the sheet 
becomes thinner. 

Introduce the following dimensionless variables: 


W = -2L- z = -Z-* x - 2L- t = w = 22 L 
W 0 ’ W 0 ’ T 0 Wo’ w 0 

Rewriting equation (2-10) in terms of the new dimensionless variables yields: 


( 2 - 11 ) 


* 




di = -. 


4x 2 


dt w o , w§ l x gW 0 t A 
W 0 g W 0 g y w 2 ) 

Or, in terms of the Froude number, defined as the ratio of inertial forces to gravitational 
forces: 


gf, One obtains: 

' -Vi 

* 


Fr = ^°_ 
gW 0 


4x J 


dt ft(l + t(£+l| 


^ For the special case of small t /Fr it is permissible to neglect it, giving: 


( 2 - 12 ) 


dl = . 


4? 


dt Fr(l-t-x) 

Rearranging and integrating results in an equation for t as a function of t. 


(2-13) 




$ 


& Inx - 1 1 + c 

tv x Fr 

■'ii Using the boundary condition that at t = 0; x = 1 it can be seen that c = -1 and: 
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(2-14) 


Writing equation (2-4) in terms of the previous dimensionless variables results in: 


tV /Cv 


I 


dW _ / 8 ct _ 

dt V wgp-Tot 

Incorporating the Weber number, defined as the ratio of inertial forces to surface tension 


sti& forcc5. 


f& one obtains. 


* ,T 


. V*52 

$ * f 






We = pWSTo 

CT 

(2-15) 

dw_ /_a_ 

dt V Wex 

(2-16) 


If equation (2-13) is divided by equation (2-16) an expression for dt/dW can be 
obtained. 


m 

^'*33 

r if 


dx dt -4-r 

dldw /X Fr (i + ^) 

V Wet 


dj_ _ m7Ez 5rl 
dW Fr ( 1 + x) 


This can be integrated to find W(x) 




IH 


3*i=m£[ d w 

5/2 Fr / 
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*t7* 


,2^-3/2. 2 --l/2 = 3 2^W + C 
3 Fr 


W = -B- (c - J x 0/2 - 2t ' 1/2 1 

VZW?l 3 / 


Using the boundary condition that at W = 1; x = 1, it is found that: 
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Therefore: 



C =£ W + £ 

Fr 3 


w = -£i-fiZWE + 2..2z-3,2 2 --m\ 
V2We v Fr 3 3 J 


W=, -3 C ^' M + S -' /2 -“) 


(2-17) 



yields: 


Writing equation (2-6) in terms of the dimensionless variables introduced earlier, 


or 


Integration yields z(t): 




dr 


W S 


f = j£-+l 

dt Fr 




dt 




t -2 - 

Z = r^t Z + t + C 
2Fr 


Applying the boundary condition, at t = 0; z = 0, results in c = 0, therefore: 




( 2 - 12 ) 


^ In order to find the sheet length over slit width ratio, IVWq, equation (2- 17) should 

H t* SOlVCd for T nume rically when W = 0. This value for ? can then be plugged into 
M, ^ ^° n ^ CO *' Finall y » f ^ e vaJue (or t can be plugged into equation (2-18) to 
'& find Z ’ 1116 rCSulting vaIue for 5 is the sheet length over slit width ratio for the given flow 


■ffs 

pi parameters. 

; lk, 


The value for IVW 0 can be determined much more easily and without the aid of the 

& PUtCr ** giaVlty is ne 2 Iected - 35 is the case in space application. As was shown before 

% 
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if g = 0 then dx/dt = 0, and consequently x = 1 throughout the sheet. Thereby equadon (2- 
16) can be written as: 


& 




dw = or 

dt V We 

This can be integrated to find W as a funcdon of t: 




dt 




W =-V4 J+C 

£ From the boundary condition, t = 0; W = 1, it can be seen that c = 1 and: 


W =-V¥‘ +I 

Now, W can be set equal to zero to solve for t, giving: 




Since, when g = 0 the Froude number is equal to infinity, equation (2-18) may be rewritten 


as: 


Therefore: 


z = t 






Therefore, the L/W 0 ratio in zero gravity, is simply equal to V 8 . 


(2-19) 


2.1.2 - Computer Solution 

The Fortran program FLUDYN has been written to compute the L/Wq ratios for a 
typical range of Froude and Weber numbers. This program computes the L/Wq ratios both 
including the effect of gravity and neglecting it. The results of this program have shown a 
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negligible difference (<0.1%) in the calculated L/W 0 ratios if gravity effects arc considered. 
For this reason, the rest of the calculation for the cylinder cross-secdonal shape shall be 
done by neglecting the force of gravity. 

The text of FLUDYN can be found in Appendix B. Figure 2-3 shows how W and 
x vary with z for a Weber number of 20. This shows how the sheet becomes thinner on the 
way down and how it keeps a basically triangular shape under the influence of gravity. 

The choice of Weber number equal to 20 was to show the sheet thinning, as the Weber 
number gets higher this effect becomes less and less. Figure 2-4 shows the sheet length 
over slit width ratios versus the Weber number. These were calculated without the effect of 
gravity. Those values that were calculated with the effect of gravity could have been 
plotted but the difference would not have been noticeable. 



Figure 2-3 Sheet Width and Thickness vs. Vertical Position for We =20.0 


3.0 


j= 


to 

5 ? n — 

/ 






j 2 - u j 

S i o - 

( 






g 1 . U i 

00 J 

0.0 A 

i 

— i — r 7 T-r 

i | * 1 1 

i i 1 i 1 1 

— r-T — T — ! — r 

1 1 l l j 

■■ » i ■■ 'i 


0.0 100.0 200.0 300.0 400.0 500.0 600.0 

Weber Number, We 

Figure 2-4 Sheet Length/Sheet Width Ratio vs. Weber Number 
2.2 - Derivation of the Edge Cylinder Cross-Sectional Shape 


The full equations of continuity and motion, including the surface tension pressure, 
arc as follows (see Appendix A): 


Continuity: 


| = -(Vpv) 

(2-2C) 

Motion: 

a 

^ V = 

- V-pvv - V(P a + Pst) - + Pg 

(2-21) 


Several assumptions need to be made to simplify these equations. 

1) The flow is steady state; 

2) The fluid is incompressible; 

3) Gravity is neglected; 



v-' < \V-„ 


4) The fluid is inviscid; 

5) Air pressure gradient is neglected; 

6) The flow is irrotadonaJ; 


As a consequence of these assumptions, 

dp d 


aT = 5 pv = 0 

p = constant 
g = x = 0 

VP a = 0 
V x v - 0 


and the governing equations reduce to, 
Continuitv: 


V- v = 0 


Motion: 


V.vv + ^I = 0 


Through the use of vector identities, the equation of continuity, ar\d the irrotationality 
assumption, the first term in the equation of motion can be simplified. 

V- vv = v- Vv + v(v- v) = v- Vv =J-V(v-v) - v x (v x v) = lv(v • v) = ^-Vv 2 


Continuity: 


V-v =0 


Motion: 


2 P 


2.2.1 - Integral Formulation of the Continuity Equation 


Consider Figure 2-5 showing a typical cross-sectional shape and the coordinate 
system to be used. The function s(x,z) describes the cross-sectional shape. The function 
r(z) is essentially the width of the edge shape, and is defined as the point where s = Tq/2 
and u = -u^. 



Figure 2-5 Typical Edge Cylinder Cross-Sectional Shape 


Expanding the continuity equation for this control volume yields: 


Jx ay 


Integrating equation (2-24) from 0 to s over y gives: 


/>!]»■” 


(2-24) 










jf lr (udy,+ / o Ir^)' 0 


The second integral in equation (2-25) may be directly integrated. Knowing that v at y = 0 
is equal to 0 due to the symmetry of the sheet, 


• (vdv) = Vj.- Vy = o = v s 


The first integral in equation (2-25) may be evaluated through the use of Leibnitz’s rule. 


1 Jr (udy)= Jr/ 0 ud y-(us^- u y =o^(°))=^s(i J o udyj-u 
If (u) is defined as the average x-component of velocity over y, then: 


-if 

Jo 


( 2 - 27 ) 


|( t;d y) = |- S ( U )- Ul | 


(2-28) 


Substituting the integrals in equations (2-26) and (2-28) into the integrated continuity 


equation, equation (2-25), gives: 


3 * i 3s a 


Since s = s(x,z) it may be written that: 


ds _ ds /dx) ■ ds /dz) 
dt 3x\dt/ 3zldt/ 


It can be seen that for this case: 


= Wj = w 0 


(2-29) 




Therefore: 


dt ± 



(2-30) 


Substituting this relation for v s into equsun^S gives: 


3 / 1 3s — - ► » 

tt— s(u) - u 4 + a* - — - - — - 0 
dx dx ' 


on 


3x 


s(u) = f- 


(2-31) 


If equation (2-3 1) is integrated from 0 to x it is f or-— 


f 


3x 


s(u)dx = -«? — x 


(2-32) 


The first integral in equation (2-32) may be dir 
(u) = s = 0. Therefore: 


, knowing that at x = 0; 


jf ^ s(ujdx = m ' 


(2-33) 


The second integral in equation (2-32) may be e var.T atuu gh the use of Liebnitz's rule. 


I sdx -l s '$- i ~P ) hdtl sdx 


If equations (2-33) and (2-34) are substituted back tKaation (2-32) it yields: 


( s ( u ))x » -wo -4- 


r- 

I 

A 







(2-34) 


(2-35) 




j Lilii 



If A is defined as the cross-sectional area of the sheet from x = 0 to x = x for any z 


then: 


= 2 f sdx 
J o 


and 


If A,, is equal to the cross-sectional area of the edge cylinder then: 


r 

Ac = 2 1 sdx 

Jo 

If the equation (2-37) is applied with x = r, refer to Figure 2-5, then: 






A = Ac 


Therefore: 


To /~2o~ _ w 0 dAc 
2 V pt 0 2 dz 


Xo^f—^—dz - dAc 

V PWqT 0 

t 0 VwF dI ‘ dA<! 

Equation (2-39) can be integrated giving: 


T °Vlr/ dz= / 


dAc 


(2-36) 

(2-37) 


(2-3S) 


(2-39) 


Ac = To Vw£ Z + C 
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^Knowing that when z = 0; = 0 it can be seen that c = 0 and: 


We 


(2-40) 


122.2 - Formulation of the Motion Equation 


The motion equation can be written for any two points on the surface of the sheet. 
jJLet the first point be any arbitrary point, s, on the surface. Let the second point be the 
•HP Stagnation point in the center of the sheet on the surface, where u = v = 0. Writing 

* JK*. ' & 


rigequation (2-23) for these two points vields, (see Appendix A): 

r + 1 
ax2 + 0 z 2 1 


+ V? + wg] - 0 


l + /^)%/^i 2 ]2 

lax/ [dzj 1 


2 P[ u stag + ''stag "*■ w o] ' O 


a^s + 3^S 1 

5x2 dz 2 


\dxl [dzj J J stag 


f a^s + 5 2 s 

5x 2 5z 2 

— . rT 

a 2 s 5 2 s 

3x2 + ^2 

U lax/ + [dzj j J 

— - u 

l 

h&'Srfj 


2 P [u| + v|] - a 


,su « (2-41) 

Since all of the partial derivatives on the right hand side of equation (2-41) are taken at the 

ir 

^Stagnation point, all of them are constant Therefore the right hand side is equal to a 
inconstant, or: 



2 pf u s + v?]-<7 


a 2 s ( 5 2 s 
5x2 + 0 Z 2 
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= —£L 
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stag 


(2-42) 
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9*^1 S^tj 

K 2 + ac 2 




1 + — 

Ui 




To 

2Rer-j 
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i: . 
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n.*7T' 


ii ^ +f H|u?r> + f-f]-5^-][ > + W* 

Hr- a? 2 a? 2 2 [ L WJ R >“4 UdJ (2 . 

From the continuity equation it was shown that: 

™jl (2 , 

. , ; ||ln order to simplify the expression for U s , it will be assumed that u s is equal to (u). This 
f - removes any y-dependence on u, giving: 


i 

* <v 


(2-43) 

-# . 
'and: 



Jr 


(2-49) 


‘ Ip; ? ^Therefore, if equations (2-48) and (2-49) are substituted into the definition of U s , given in 


^equation (2-46), it can be seen that: 

■ i> 

-i r 


fdA) MA) 

_ tq Idzq _ l ' dz/ x 
2s dAc q dAc 
dz dz 


(2-50) 


• If it is assumed that the cross-sectional shape keeps a similar form for all z, then the 
definitions for the cross-sectional area of the sheet and the cross-sectional area of the edge 

V* 0 

flinders may be scaled by r, on 

■ t 
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ftf 2.2.3 - Three Possible Simplified Solutions 

The problem still remains to solve fort) as a function of \ and £. It would be 




mi 


& 


m 


m 

T easiest if all the z-dependence were removed. That is, to be able to solve for q as a 
% function of % for any given A c . The only place where any z-dependence remains in this 
equation is in the term. Three simplified solutions to this problem will be 

considered. In the first simplified solution, (i), this term will be neglected. In the second 
simplified solution, (ii), this term will be set equal to a constant. In the third simplified 
solution, (iii), this term is set proportional to the other second partial derivative term, that 
is, 

0^-0 

ac 

a^T\ 

ii) = Q = constant 

oC 2 

..., d 2 !} 

m) — - = a — L 

3C 2 ac 

Since, cases (i) and (ii) can be regarded as limiting cases of the' more general case (iii), the 
solution for this case will be considered first 

Taking assumption (iii) and substituting into equation (2-58) yields: 


£i + a£a-i 

d? 2 d$ 2 2 




ifl 

timAcM U/J 


to 


‘‘Stag 


Wl 


f 


(2-59) 


This is a second order, non-linear, ordinary, differential equation which may be solved for 
q as a function of £ for any A c . The boundary conditions on this equation arc as follows: 


*1 


1) at 2; = 0, q = 0 and — - — ><*» 

dC 

2) at q = r, q = 1 and — = 0 

dC 


The solution can be found by using a fourth order Runge-Kutta technique. The 
algorithm uses what can be thought of as a shooting technique. and a are free 

variables. Knowing A c , XQ/R SU g and a are selected. Beginning at q = 0, where T| = 0, the 
algorithm calculates T| at At; and from that calculates rj at 2Aq and so on. The program is 
allowed to shoot until r| passes through a minimum. This point is f. This insures a zero 
slope where the edge cylinder connects to the sheet. Thus the first part of boundary 
condition 1 and the second part of boundary condition 2 will always be satisfied. The next 
step in the routine is to insure that the proper tf/R^g sad cc were selected. 

First, the value of r\ at q = r is checked to see if it is equal to 1. If it is not, the 
value for a is changed and the algorithm is started again. This continues until the value of 
•q at r is equal to 1. Second, the value of A at q = F is checked to see if it is equal to A c . If 
it is not, the value for x^/R^g is changed and the algorithm is started over again. This 
most often results in a being changed again, but eventually both conditions are satisfied. 
This satisfies the first part of boundary condition 2 and boundary condition 3. 

Due to the second part of boundary condition I, the second point, r|(A£), must be 
found by an analytic approximation. For small A is approximately equal to zero, and 
equation (2-59) reduces to: 


d 2 !) 

d § 2 


xo 


1 + 


(dn) 2 ! 

k/J 

2Rstag ( 1 "*■ ®) 


(2-60) 


The solution to this equation is: 



30 


■■l'XrA'7 



[ 2 ( l+a)R JU , e 

dri 1 Tn * S 


(2-62) 


The Fortran program SHAPE has been written to solve equation (2-59) using these 
boundary conditions and the approximations for small 2;. The listing of this program can 
be found in Appendix B. The results of running the SHAPE program for a number of 
Ac’s are shown graphically in Figure 2-6. Figure 2-7 shows that there are only specific 
values of a and that can exist for any given A c . Tnese two quantities quickly 

diverge as A c approaches 0. For this reason, results for A c 's less then 40 were unable to 
be obtained- One interesting note is that if assumption (i) had been used it would have been 
equivalent to setting a equal to 0 in this solution. As can be seen from this graph, the only 
A c at which the boundary conditions could be satisfied had this term been neglected is 
approximately 454. For this reason it is apparent that assumption (i) would have been 
invalid- Figure 2-8 shows how the quantities Tima* and ? vary with A c . T),^ is the largest 
value T) reaches in the cross-sectional shape and represents half of the maximum thickness 
of the cylinder. Figure 2-9 shows the ellipticity of the end cylinder. This is particularly 
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Figure 2-6 Cross-Sectional Shape of the Edge Cylinder for Case (iii) 
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! Figure 2-9 Elliptirity vs. Dimensionless Cross-Sectional Area for Case (iii) 
^in te r esting in that it shows for small areas the cylinder is wider then it is thick, but for 

m 

^larger areas the cylinder approaches a more circular cross-section, which would have been 
i'Aubtintively expected 

.i^-1 Taking assumption (ii) and substituting into equation (2-58) yields: 


d^ 2 2 [ti 2 Ia c /[ 


1+ — Z -Q 

U/l 


(2-63 


te-Tnis equation is solved similarly to the previous one. The method is the same and the 
S&v. 

^boundary conditions are the same, the only difference is the two free constants are now 
^Jo^stag Q- The second point ti(A^) still has to be found using an analytic 
^approximation. If A and Q are set equal to zero in equation (2-63) it yields: 
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Figure 2-11 T^/Rc^g and Q vs. Dimensionless Cross-Sectional Area for Case (ii) 


.S'- 


0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 3500.0 

Dimensionless Cross-Sectional Area, A 


Figure 2-12 and r vs. Dimensionless Cross-Sectional Area for Case (ii) 
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same boundary conditions as before and these approximations for small Jj. The listing of 
this program is also appended. Figure 2-10 shows the results of SHAPE2 when it was ran 
for a number of A c 's. Figure 2-11 shows there is only one specific tc/R^tag and Q for any 
given A c . This is also in agreement with the previous solution in that if assumption (i) was 
used the only value of A c where the boundary conditions could have been satisfied was 
454. Figure 2-12 shows how t] max and r vary.\yith A c . This is quite peculiar where it 
shows that for A c ’s less then 497 a decrease in area actually produces an increase in F. This 
is definitely not to be expected. The computer solution became quite unstable for small 
values of A c , for this reason, no values were obtained for A c ’s less then 400. Figure 2-13 
shows the ellipdcity using assumption (ii). This is quite unexpected, for A c ’s greater 



Figure 2-13 Ellipticity vs. Dimensionless Cross-Sectional Area for Case (ii) 

then 497 the cylinder becomes thicker than it is wide. This is very apparent in Figure 2-10 
where the curve corresponding to an A c of 497 is dashed and the curves below it have a 










large F. This is very different then case (iii) where the cylinder approached a circular shape 
; - for large A c 's. For these reasons assumption (ii) will have to be discarded and assumption 
(iii) accepted as the best of the three simplified solutions. 

With assumption (iii) being accepted as a physically possible solution, all that 
remains is to determine A c and run the SHAPE program for this value. It was shown 
earlier than 


?§ but ’ 

- ' 

•^.Therefore: 

h. 


Ac = ^-Ac 
*0 


Ac = tozy 


HI 

We 




(2-56) 


(2-40) 


(2-67) 


■from Chapter 2, Section 1 it was shown that: 

W*- 


if 

:$L 


Wo V 8 


(2-19) 


•Substituting equation (2-19) into equation (2-67) gives: 

rrT* * 






XqL 


(2-63) 


^With this expression, all that has to be known is the dimensions of the slit, the length of the 
jjfMt, and how far down the sheet the shape is desired to be known to calculate A c . Once 
*.lhis is known, the program SHAPE is ran and the output describes the theoretical cross- 
S ecti onal shape. 



Chapter 3 
Stability Analysis 


The focus of this chapter is to analytically predict the stability of a non-planer, thin 
liquid sheet. It has already been shown by Geoffrey Taylor [1] that a perfectly planer 
liquid sheet is completely stable. It has also been shown by Rayleigh [4] that a cylindrical 
jet is unstable to an asymmetric disturbance. Unlike Taylor’s work, the sheet in this study 
is not perfectly planer and does not have a zero velocity throughout in the transverse 
direction. The sheet has curvature, particularly in the edge cylinders and also in the sheet 
itself due to slit defects. With this in mind, this particular analysis may be considered a 
combination of Taylor's and Rayleigh's work. 

By neglecting gravity the velocity field to be considered is the following: 

V =e> u(x,y,z), v(x,y,z), w 0 

Without loss of generality a Galilean transformation may be made to reduce the z-direction 
velocity to rest, leaving: 

V = u(x,y,t), v(x,y,t); w = 0 

Assuming an incompressible fluid, the continuity equation may be written for this flow 
field as: 

3u 3v . 

K + JT° 

Assuming irrotational flow, it may be written that: 


9u 9v _ Q 
9y 9x 

Introducing a velocity potential to satisfy the irrotational assumption, 
















:’ v; •jjjc flow on the sheet surface must satisfy equation (3-2) and the boundary conditions, 
‘ l equao c,ns (3-3) and (3-4). 



If <£o and Sq are the solutions to these equations, a small perturbation or disturbance 
and As may be added to these solutions as follows: 

<f> = $o(x,y) + A<J>(x,y,t) ( 3 . 5 ) 

s = so(x) + As(x,t) ( 3 - 6 ) 

: values for o and s may be substituted into equation ( 3 - 2 ) to obtain the disturbance 
equation. 



— - (do + Ao) + ^ — (<?0 + Ad) = 0 
dx z dy z 


d 2 pp d 2 Ao d 2 do d 2 A6 

dx 2 dx- dy 2 dy 2 

jtemoving the terms that are identically zero by equation (3-2) leaves: 




tZ: 


d 2 A6 d 2 A<J> 

- + = 0 

dx 2 dy 2 


(3-7) 


■ ijg Equations (3-5) and (3-6) may also be substituted into the boundary conditions, 
i^pbstituting equations (3-5) and (3-6) into equation (3-3) yields the following. 


4 

or. 


Jp- (<t >0 + A<j)) s = Jj- (so + As) + ^ (sp + As) ^ (<|)o + A4>) s 


^dA<|)' 

\3y It [dy J s ~ dt ' dt 


dso + dAs 


dsp f d0o ) dAs f d<}> 0 ) | dsp f dA(ft | dAs (dA<b \ 

dx [dx J s + dx \ dx ) s dx [ dx ) s dx [ dx / s _ 


terms that are identically zero by equation (3-3) can be removed from the above 


,.^r 

L#.. 


By the assumption that A<j) and As are small, the last term in the above equation 

■c 

ff^y be neglected, thus resulting in a linear equation in A<1> and As. 
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• •»*■ ••• . ' • ■•/ ..; -4 -V.r. 


• “0 

As = sefc* • “0 


where to is the frequency and k is given by the following relation: 


k = 2i 

X 


(3-10) 

(3-1!) 


(3-12) 


where X is the wavelength. These wavelike forms may be substituted into the disturbance 
equation, equation (3-7), and into the boundary conditions, equations (3-8) and (3-9). The 
resulting equations may then be solved for co. If there is a positive imaginary part to co then 
any initially s m a ll disturbance will grow exponentially with time causing the solution to be 


unstable. 


These wavelike forms can be substituted into equation (3-7), yielding: 




JL. (<f>ei(kx - cot)] + - cot)) _ o 

ax 2 dy 2 
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or. 


(^ikpe^ - ®0 + e^ - “ l )— - = 0 , 


£-*-o 

dy 2 


(3-13) 


4he solution to equation (3-13) is as follows: 

W 


<{> = cisinh(ky) + C2C0sh(ky) 

V*‘ 

where cj and C 2 are constants. Therefore, substituting equations (3-14) and (3-10) 


(3-14) 


..Ration (3-5) gives: 

.. •0' 


<{> = <f>0 + (cisinh(ky) + C2Cosh(ky))e i ( lc * - “0 


(3-15) 


Since every wavelike disturbance can be considered as a combination of a 

WnwY 

’Proetric and an asymmetric wave, these are the only two waves that need to be 
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considered. In the symmetric case, the y-velocity on the back surface of the sheet is the 
opposite of that on the front surface. In the asymmetric case, the v-veiocity on the back 
surface of the sheet is equal to that on the front. These two cases are shown in Figure 3-1. 
The boundary conditions for these two cases are as follows. 

1) Symmetric Case 

v s = '-Y-s 


ctyj j d<b 

dyl~ wy 


( 3 - 16 ) 



Symmetric Case 



Asymmetric Case 

Figure 3-1 Schematic of the Symmetric and Asymmetric Waves 


2) Asymmetric Case 


v s = v -s 


: -Mr' 


43 







wy/s Uy L (3-n; 

For case 1, the symmetric case, boundary condition (3-16) can be applied to 
equation (3-15), giving: 

(liy^) + ( ClCOS ^( ks ) + C2sinh(ks))e i ( ,c ' -«*) = - * (ciCosh(-ks) + C2Sinh(-ks))e l ( kx * ox ) 

The first terms on both sides of this equation are equal from equation (3-16). Rewriting 
using trigonometric identities gives: 



(cicosh(ks) + C 2 sinh(ks))e^ la ' 0M ) = - (cicosh(ks) - C2sinn(ks))e i ( kx • mt ) 


Simplifying yields: 


c 1 = -ci 


Therefore, for the symmetric case, Cj must be equal to zero. Now if c 2 is set equal to A^, 
the amplitude can be written as: 


<J) = A S yCOsh(ky) , Symmetric Case (3-18) 
For case 2, the asymmetric case, boundary condition (3-17) can be applied to 
equation (3-15), yielding: 

(ip) + ( ClCOsh ( ks ) + C 2 sinh(ks))e i ( kx - “*) = + (cicosh(-ks) + C2sinh(-ks))e ! ( kx * “0 

According to boundary condition (3-17) the first term on both sides of this equation are 
equal. Rewriting with the use of trigonometric identities gives: 


(cicosh(ks) + C 2 sinh(ks))e i ( l01 - “0 = (cjcosh(ks) - C2sinh(ks))e i 0« * °*) 


Simplifying gives: 


c 2 = ‘ c 2 
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Therefore, for the asymmetric case, Ct must be equal to zero and if Cj is set equal to A^, it 
yields: 

^ = AasySinhCky) Asymmetric Case (3-19) 
The wavelike forms may also be applied to the boundary conditions. If equations 
(3-10) and (3-11) are applied to equation (3-8) it gives: 

e i(!ct-Qx) |i| =se^- mt \-ico)+ se‘( ls ^' <nt \ik)|^j +^^ s c^' CM \)k) 

Dividing out the exponentials and multiplying by i gives: 


[3d I *Jd<bo\ 9so?, 

la7r s “- s W-a^ 


(3-20) 


If equations (3-10) and (3-1 1) are applied to equation (3-9) it yields: 


*■' w + |a y ^a y / i e P( 1+ /awfp /2 

Ux/ J 


Simplifying yields: 
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(3-21) 


Simultaneously substituting equation (3-18), for the symmetric case, into both 
boundary conditions, equations (3-20) and (3-21), at y equal to Sq yields: 


ikA sy sinh(kso) = sco - u 5 ,osk - |^jkA sy cosh(k;so) 


(3-22) 


ana 

- icoA^yCOsl^kso) + iku,,oA,yCosh(kso) + kv 5 ,oA 5 ySinh(kso) +^- — rr 

\ \ dxl) (3.23) 

Solving equation (3-22) for s and substituting into equation (3-23) yields: 
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co' = to - ku^o = co 




Therefore, equations (3-26) and (3-29) reduce to: 


co“ - ~~ tanh|^-j = 0 


Symmetric Case (3-30) 


k*^CT 

or - cot 




Asymmetric Case (3-31) 


These two equations clearly have two real roots each. Therefore, this analysis predicts 
complete stability for any disturbance in the flat area of the sheet This is in agreement with 
the work done by Taylor. It is possible in this case to solve these equations for co and 
consequently predict the phase velocity for each case. 

Solving equations (3-30) and (3-31) for co yields: 


co = -y ^2. tanhj^j 


Symmetric Case (3-32) 


=V1M¥) 


Thus, the wave velocities can be calculated as: 


Csy = k = Vo ^2) 


c 

c “ y_ k - 


^-coth|“ 


Asymmetric Case (3-33) 


where c sy and c^ are the symmetric and asymmetric wave velocities, respectively. Taking 
the limit as x goes to zero results in: 
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clcngth- The velocity of 


The velocity of the symmetrical waves is still a function of the wav 

the asymmetrical waves is independent of the wavelength and equal t° di 

1 In *i shccc cno^n^ 

Thus the asymmetrical waves are propagated wfthout any dispersion* 

f constant phase are at 

with velocity w 0 the asymmetric waves will appear at rest if tines 01 
an angle V * to the direction of flow, provided: 


sinV'-^ 

But if the angle the sheet edge makes with the z-direction is V then. 


(3-36) 


ti- 
tan w = — 

Y wo 

Thus, for a long sheet, high wq, these two relations can be simplified to* 


Y W 0 


and 
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wq 


(3-37) 


, edge. The 

Therefore, for a long sheet the asymmetric waves are parallel to ** 

f j disturbance. 

velocity of the symmetric waves is still a function of the wavcleng^' ** 

Therefore the symmetrical waves appear as a set of parabolas on tfv? 

i* figure 3-2 waves 

Figures 3-2 and 3-3 show actual photographs of the sheet, /" ^ . tf he 

clearly visible that are parallel to the edges of the sheet. These are feff***** 13 * 
asymmetric waves. In Figure 3-3, there is clearly a set of parabolic 

waves, 

emanating from a single point. These are representative of the symr,^ 7 ' 







Substituting equation (3-38) into equations (3-26) and (3-29) results in: 


(co'}~ + iku 5< o^tanh(kso)co' + ( i^. . tanh(kso)] = 0 

3X P' l3x 1 SymmenicO* S3, 


(o') 2 + ikUs.o^cothlksoJoj' + ^fi^-- coch(kso)} = 0 

P X _ Asymmetric Case C-C 1 .- 

Rewriting equations (3-39) and (3-40) in terms of the following variables 


X = kso.y = ^.Y=— ? 


Us, ° ‘ pHqSo 


results in: 


Y 2 + i^Xtanh(X)Y + yX 3 (i^. - tanh(X)l = 0 

™ l ox I Symmetric Case G-O 


Y 2 + i^Xcoth(X)Y + yX 3 (il^- - coth(X)j = 0 


Asymmetric Case (3-43) 
Taking the limit of equation (3-42) as X, the thickness, goes to zero results hr 



or, 


Y 2 = 0 Symmetric Case (3-44) 


Taking the limit of equation (3-43) as X goes to zero results in: 



or. 


Y = - 


. 3sp 

*8x 


Asymmetric Case (3-45) 





53 


Chapter 4 

Description of Experimental Apparatus and Study of Hole Formation 

A test rig for producing a liquid sheet of silicon oil in a vacuum at NASA Lewis 
Research Center was constructed previous to this study. This facility was constructed to 
study the fluid dynamics and emissivity of a chin liquid sheet. This test rig produces a 
relatively large liquid sheet, on the order of ten feet long. The sheet was contained in a 
solid steel tank in order to hold a vacuum. The only way of viewing the sheet was to either 
look through small windows at the top of the tank or look at the videos produced by a small 
camera with the ability to move up and down and around inside the tank, either way only a 
small portion of the sheet was visible at any one time. 

The sheets formed in this rig appeared to be unstable. That is, when the sheet was 
observed it would appear to flicker. The sheet also would splatter near the bottom, making 
it virtually impossible to place a probe near the bottom of the sheet to take any temperature 
measurements to determine the emissivity. An investigation was undertaken in order to 
determine the source of this instability. High speed film loops were taken of the sheet. 
Observation of these high speed film loops revealed that holes formed in the sheet and that 
they grew as they went down the sheet. 

In order to further investigate the formation of these holes and to study liquid sheet 
flows, three small scale test rigs were made. The sheets in these test rigs were formed with 
water, in air, and under the influence of gravity. The main advantage the these test rigs 
was that the entire sheet was visible at all times. 


4.1 - Operation of the Small Scale Test Rigs 


A schematic of the general design of the small scale test rigs is shown in Figure 4- 
1. The test rigs basically consist of a plenum or tank with a small slit in the bottom. When 
the tank is pressurized with water, the water is forced through the slit and the sheet forms. 


Air Escape and Water Overflow 


Water In 



Figure 4-1 Schematic of the Small Scale Test Rig 
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The water from the sheet is then collected in a bucket or water collector which is equipped 
with a drain. 

Tap water flows into the system through Valve A. Valve A is used to control the 
mass flow rate through the system and the pressure in the plenum. After the water flows 
through Valve A, the water is ran through a filter rated to remove any particulates larger 
than 10 jam. During the time the plenum is filling up with water. Valve C is shut and 
Valves B and D are open. Valve B allows for air to escape from the plenum during the fill- 
up time. Once the tank is full, water flows out through Valve B. At this time, Valve B is 
closed and the plenum begins to pressurize and a liquid sheet begins to form. The water 
from the overflow and the liquid sheet is collected in the water collector. Valve D is almost 
always left open, when it is closed water collects in the bucket and mass flow rate 
measurements can be obtained. When the testing is over, Valve A is closed and Valves B 
and C are opened to allow for w'ater to drain out of the tank and air to enter the tank. 

4.2 - Differences in the Three Small Scale Test Rias 


The differences in the three test rigs are due to the size of the plenums and the size 
of the slits that they run with. In the first test rig, the volume of the plenum was 66 in^. 
The slits all had widths of 1.345 in. and thicknesses of either 0.0021 in., 0.0032 in., or 
0.0040 in. The second test rig used the same three slits as the first, but the volume of its 
plenum was 1,900 in. The third test rig had the unique feature of an adjustable thickness 
on the slit. The slits on this test rig were 8 in. long and the thicknesses were adjustable 
from 0.0200 in. to 0.0015 in. Also, an impediment could be placed in the slit in order to 
form sheets less than 8 in. wide. The volume on this plenum was 3,910 in^. 

Figure 4-2 is a photograph of the third small scale test rig, showing the large 
plenum and the sheet at the base of the plenum. Figure 4-3 shows a close-up photograph 
of the adjustable slit from the third test rig. The slit is 8 in. long and the positioner has a 




Figure 4-2 Photograph of the Third Small Scale Test Rig 


digital readout, capable of showing the thickness of the slit to accuracy of 0.00005 in. 
Since it was necessary to produce both sides of the slit 8 in. long, parallel, and :o the 
highest possible degree of precision there was a great amount of time spent in the 


machining of this slit. 
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Figure 4-3 Photograph of the Adjustable Slit 


4.3 - Experimental Procedure for Determining the Surface Area Lost Due 
to Hole Formation 



The first test rig, equipped with the 0.0032 in. thick slit, was used for determining 
the surface area lost due to hole formation. The photographic technique employed for 
determining the surface area lost due to hole formation is shown in Figure 4-4. A camera 
was set up directly in front of the sheet and a strobe light was set up off center behind the 
sheet, the dashed line represents the line of sight of the camera. The strobe was placed off 
center so the area directly behind the sheet would not be flooded with light. The camera 



used was a Pentax KX and the film was Kodak TMAX400. The shutter speed on the 
camera was set to 1/60 of a second and the strobe was set to 18,000 flashes per minute. 
This resulted in there being 5 images of the sheet on each frame of film. Therefore, i: v/as 
possible to see an image of a hole near the top of the sheet and see four consecutive images 



v , 

yi# of the same hole further down the sheet all on one photograph. Figure 4-5 shows a tvpical 



Strobe Light 


Za 



Camera 

Figure 4-4 Top View of the Experimental Setup for Hole Formation 

r 

one of these photographs, showing an image of a hole near the top of the sheet and three 
more consecutive images of the same hole as it goes through the sheet 
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Figure 4-5 Photograph of a Sheet Experiencing Hole Formation 
Slides 7 / ere made of each of these photographs. These slides were then projected 
on to a digitizing pad. The outlines of the sheet and the holes were then digitized and Lie 
areas of the sheet and the holes were computed using the Autocade software package. 
From these areas the percent of surface area lost could be calculated. This was done for a 
number of sheet lengths and operating pressures, 'flic results are shown in Figure 4-6. 


Figure 4-6 clearly shows the random appearance of holes in the sheet. Even though 
there was a general increase in surface area lost as the pressure was increased, there also 
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Figure 4-6 Surface Area Lost vs. Line Pressure for the First Test Rig, T = C.0032 in. 
were distinct exceptions to this rule. This led to the hypothesis that the cause of the holes 
was not directly linked to the fluid dynamics of the sheet flow alone. In addition, the holes 
never exhibited any pattern for where they appeared. Holes appeared on the left, right, and 
center of the sheet. Due to the apparent random nature of their occurrence, it was suspected 
that upstream disturbances might be the cause of the holes. 

For this reason, the second test rig was constructed. With such a large plenum it 
was assured that the water above the slit would be nearly stagnant and the pressure 
fluctuations would be greatly minimized. When photographs were taken of the sheet with 
the new plenum anached it was revealed that no holes were ever formed for all the 



conditions tested. 


Thus it was concluded that upstream fluctuations were the cause of the hole 
formation. This finding necessitated the construction of a new plenum for the test rig at 
NASA Lewis. Since the test rig at NASA was undergoing renovations the experiments 
performed on the liquid sheet and the data presented in the following chapters were 
generated by the second and third small scale test rigs. 


n / 




Chapter 5 

Experimental Investigation of Liquid Sheet Flow 


Three interesting experimental studies were conducted on liquid sheet flows. In ti 
first, the sheet length over slit width ratio was determined for a number of slit sizes and e: 
velocities. The results from this could be compared to length over width values calcuintec 
in Chapter 2, Section 1. In the second study, the widths and thicknesses of the edge 
cylinders were measured and compared to the computer generated results of Chapter 2, 
Section 2. In the third study the effects of air resistance on the sheet was studied. 


V- 
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5.1 - Determining the Sheet Length Over Slit Width Ratio 

The second test rig, equipped with the 0.0021 in. and 0.0032 in. slits, was used fr 
calculating the sheet length over slit width ratio; Since the water collector on the test rig 
was set up with a valve on the drain the only equipment necessary for the calculation of the 
sheet length over the slit width ratio was a millimeter scale. The water collector was 
graduated with two lines. The volume of water between the two lines was equal to five 
liters. When, a steady sheet was produced. Valve D (see Figure 4-1) on the water collector 
was closed. Once the water reached the first line on the inside of the water collector a 
stopwatch was started. While the bucket was filling up, the length of the sheet was 
measured with the millimeter scale. When the height of the water reached the second line, 
the stopwatch was turned off. Knowing the amount of water collected and the time it took 
to collect it, the volumetric flow rate could be calculated. Thus, the sheet length over slit 
width ratio could be determined for various flow rates. 
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5.2 - Determining the Edge Cylinder Cross-Sectional Shape 


5.2.1 - Description of Experimental Setup 


The photographic technique employed is shown in Figure 5-1. The dashed line 
represents the line of sight of the camera. A beam splitter and two front faced mirrors were 
used in order to view both the front and side views of the sheet on the same photograph. 
The tip of a wire was placed near the sheet edge in order to insure that both images on the 
photograph were at the same point on the cylinder. Due to the sheet going away from the 
camera in the side view, there was only one place in this view that was in focus. As a 
result, the camera was set on the tripod with the ability to move forward and back. 

Realizing that if the focus on the camera is not adjusted, then any place that is in focus will 
have the same degree of magnification as any other place that is in fccus. It was then 
possible to take a picture of the front view in focus, tape the lens to insure a constant focal 
length, and move the camera forwards and back in order to get a number of points along 


the edge view of the sheet in focus. 




Figure 5-1 Top View of the Photographic Technique for the Edge Cylinders 




The camera used was a Nikon F4 equipped with a fully extended bellows and a 
200mm Nikon lens, with a +2 diapter close-up lens. Kodak TMAX 400 black and white 
film was used. The second test rig, equipped with the 0.CO21 in., 0.0032 in., and 0.004C 
in. slits, was used. One set of data was collected for each slit. 

5.2.2 - Experimental Procedure 

Once a full, steady liquid sheet was formed, the length of the sheet was measured. 
The wire dp was placed near the edge cylinder, close to the top of the sheet and the distance 
from the tip of the wire to the top of the sheet was measured. The camera and mirrors were 
placed so the wire would have the same length in each view near the bottom of the frame. 
Two pictures were taken snowing both views, and the lens was fixed in place to insure an 
equal degree of magnification for each picture. Two pictures were then taken showing just 
the front view for greater clarity. Four or five pictures were taken showing different points 
in the edge view in focus. This procedure yielded four or five data points along 
approximately one inch of the edge cylinder. 

At this point, a picrure of a millimeter scale was taken in order to determine the 
degree of magnification of the pictures after developing. Then the wire was moved down 
allowing the top of the next set of photographs to be at the same point as the bottom of the 
previous set. The height of the wire tip was measured again and the next set of pictures 
taken. A picture of the millimeter scale was taken after each set of pictures to insure the 
degree of magnification was not changing. The length of the sheet was periodically 
measured to insure that the length of the sheet was staying constant. 

Enlarged photographs were made for taking the r and T) max measurements. The 
pictures of the scale were measured with a millimeter scale to determine the degree of 
magnification. Knowing the position of the wire tip, the blown-up photograph of the scale 
could be used to determine the height of the point that appeared to be in focus in each of the 
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edge shots. The cylinder widths and thicknesses were measured to a resolution of a 64th 
of an inch. 

5.3 - Determining the Effects of Air Resistance on the Sheet 

The third test rig, with the adjustable slit, was used for these experiments. It was 
apparent that for some of the sheets formed by this test rig air resistance would have a very 
significant affect. In cases where the sheet velocity was high and the sheet thickness was 
thin, the sheet would catastrophically destruct and form droplets, see Figure 5-2. An 
attempt w as made to determine under what conditions this would happen and to what extent 
it would happen. The critical length, L^, is defined as the length from the top of the sheet 
to the point where the sheet destructs. 




Once a sheet was formed it was noted if it appeared to destruct due to air resistance 
or not. If a sheet was fully in tact, its length, L, and thickness, Tq, were recorded. If die 
sheet was not fully in tact, its critical length, Lqj, would be recorded in addition to its 
length and thickness. This was done for both 8 and 6 inch wide sheets and a great number 
of thicknesses and lengths. 
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Chapter 6 

Experimental Results 
6.1 - Sheet length Over Slit Width Ratio Results 

Knowing the volumetric flow rate and the dimensions of the slit, the exit velocity c 
the sheet, wq (assumed to be uniform), could be calculated by dividing the flow rate by the 
area of the slit. Assuming the density of water to be l.COOkg/hr and the coefficient of 
surface tension to be 0.0728 N/m, [8], the Weber number could be calculated through the 
use of equation (2-15). Knowing the dimensions of the slit and the length of the sheet, the 
sheet length over slit width ratio could be calculated. Tne results for the sheet length over 
slit width ratio are plotted against the Weber number in Figure 6-1. The solid line 
represents the analytical prediction given by equation (2- 1 9). The circular data points are 
from the 0.0021 in. slit and the square data points are from the 0.C032 in. slit. 

6.2 - Edge Cross-Sectional Shape Results 

Knowing the dimensions of the slit, the total length of the sheet, and the height at 
which the data points were taken, the dimensionless cross-sectional area of the edge 
cylinder, A c , could be calculated using equation (2-68). Knowing the degree of 
magnification and the dimensions of the slit, ri max and r could be calculated using equation 
(2-46). These data and the theoretical predictions are shown graphically in Figures 6-2 
through 6-4. 
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Wefcer Number, We 


Figure 6-1 Sheet length over Slit Width Ratio Vs. Weber Number 

Figures 6-2, 6-3, and 6-4 reveal a very interesting phenomenon. The cylinder 
appears to be oscillating. The cylinder will quickly flatten out and then reform in its 
basically circular shape. When the cylinder is in its circular shape, the shape is in good 
agreement with the theoretical predictions. The details of this cycle are shown in Figures 6- 
5 and 6-6. Figure 6-5 shows a cross-section of the cylinder as one moves down in the z- 
dircction. Figure 6-6 shows the actual photographs of the cylinder. 

The cylinder begins in its nearly circular shape, shown in Figure 6-6a. Due to the 
high curvature where the cylinder connects to the sheet, the cylinder quickly flattens cut, 
shown near the top of Figure 6-6b. This appears in Figure 6-5 as an increase in r and a 


and r, Dimensionles 
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decrease in Instead of simply regaining its circular shape, the cylinder goes throuar 
a rebuilding process where the cylinder has a peanut-like shape. As the area of the cyiintie 
increases ^lmax uicreases but on the inside of the previous cylinder. This can be seen best 
third of the way down the edge view of Figure 6-6b, where the inner part of the cylinder 
can be seen behind the outer part of the cylinder. It cannot be seen whether the cylinder 
necks down between the two parts or not. This inner part of the cylinder grows and 
engulfs the outer part to return to its circular shape. Thus appean in Figure 6-5 as an 
increase in r) max and a decrease in r. This is most noticeable in the 0.0040in. slit. 



Figure 6-2 rj max and r vs. Cross-Sectional Area 
Theoretical and Experimental for the 0.0021 in. Slit 
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Initial, Nearly Circular 
Cross-Secdonal Shape 


Flattened Out 
Cross-Secdonal Shape 



Rebuilding or Peanut Like 
Cross-Sectional Shape 


Larger, Rebuilt, Nearly Circular 
Cross-Sectional Shape 


Figure 6-5 Changes in the Cross-Sectional Shape Proceeding in the z-Directicn 











6.3 - Air Resistance Effects Results 




The results for the air resistance effects are shown in Figure 6-7. In Figure 6-7, the 
sheet length and thickness have each been scaled by the slit width. The critical length has 
been scaled by the sheet length, in essence giving the percentage of the sheet length that is 
still in tact. Any combination of W(/tq and L/Wq that is below or on the Lq^/L = 1 line 
was found to be completely in tact. Any combination above that line was found to be in 
varied degrees destroyed due to air resistance. It can clearly be seen that in order to run a 
sheet at a high width to thickness ratio in air one needs to have a shorter sheet. 



Figure 6-7 Contour Plot for L^r/L 
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Chapter 7 
Conclusion 





Theoretical and experimental studies of the fluid dynamics and stability of thin 
liquid sheet flows have been carried out in this thesis. It was shown that this was a surface 
tension driven flow. Due to the surface tension, the sheet coalesces to a point after a length 
L and, due to the conservation of mass, the water collects in the edge cylinders giving the 
sheet a basically triangular shape. It was shown that gravity had a negligible effect on the 
shape of the sheet. 

A relation was developed linking the sheet length over slit width ratio to the 
operating parameters of the sheet, namely: 




where: 


We = 


pwjfro 


Thus the sheet length was found to be proportional to the outlet velocity of the slit. Many 
experiments were carried out which demonstrated that this relation correctly represents 
what occurs in thin liquid sheet flow. 

A computer program was developed to determine the cross-sectional shape of the 
edge cylinder given the cross-sectional area of the edge cylinder. Experiments were earned 
out revealing that the edge cylinder oscillates and that the computer solution was only valid 
for the limiting case that the edge cylinder was in its basically circular shape. 

A stability analysis was performed on a non-planer liquid sheet. This analysis 


nt*. 








showed the existence of two different wave formations on the center, flat portion of the 
sheet. Photographic evidence of the existence of these two wave formations was also 
found. It was determined that the sheet may be unstable where the y-direction velocity on 
the surface is negative, see Figure 3-2. It was shown that this is a very likely condition on 
the edge cylinder near where it connects to the sheet 

A study was conducted to determine the effects of air resistance on the sheet. It 
was shown that the longer the sheet is, the thicker it needs to be in order to stay in tact. 

As to the formation of holes in the sheet, as found at the NASA Lewis test ria, it 
was found that the holes form due to fluctuations in pressure upstream from the slit. When 
a large plenum is placed above the slit, hole formation was eliminated. 


V .... 

r-j 
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Appendix A 

Derivation of the Total Pressure on a Free Curved Fluid Surface 

Including Surface Tension 

Figure A-l shows a differential fluid element on the surface of a fluid body. The 
lengths of the two sides are dSy^ and ds^. The surface is curved, and the radii of curvature 
are shown as Ry^ and R^. The coeficient of surface tension is a, therefore, the forces 
connecting this element to the stirrounding elements are crdSy X and crds^. Due to the 
curvature of the element these surface tension forces create a force normal to the surface. 
Conducting a force balance in the normal direction to the surface results in: 


odsyx 



CtiSyz 


Figure A-l Differential Surface Element 
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. d9 VI . d0yx 

PdSy X ds yz = P a ds yx dsyz - 2ads yz sin — -j— - 2adsy*sin 


Applying the small angle approximation: 


. d0 V x d9 vx 
sin =— = --- - - 
2 2 

. d9 VZ d0yz 


results in: 


I d9 V x j d9 vz 

PdSyzdSyz = P 4 dSy X dSyz * 2(7 |dSy Z ~ + CiSyX ~2~ i 


Dividing by ds^ds^ results in: 


P = P, - a 


d9y X d9yz 

dSyx dSyz 


V : % 

; :* T 


J 

I f 

* -e; 

« ;5s 


P = P*-a 


(rL + Ryz) 


Since the radii of curvature are defined as [9]: 


Ryx Ryz 


dh t 3 2 s 
9x 2 9z 2 

, IdsV Id sff /2 
1 + 5x 9zJ ) 


die total pressure is equal to: 


P = P a - a 


3 2 s [ 3 2 s 
9x 2 9z 2 

, Ids ] 2 Ids ] 2 ?' 2 
1 + 9x (dz ) 


Define the surface tension pressure as: 


PST = ' CT 



8 2 s d 2 s 
9x 2 9z 2 

(ds) 2 Id s\ 2 \ 3/2 
[3x/ + l9zj / 
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Program Listin 
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I; m 

l: f Ml 


m 
* 1 

V ' V . 

3 : :«* 


JSf : 

# -, v - 


.,> t 

rW- '^r 


**********,.,*»*»,, *************♦*******»*******•*»********»******. 

* AUTHOR: MATTHEW S. MCMASTER 

* DATE: 12-5-91 

* 

* TITLE: FLUDYN 

* 

* PURPOSE: DETERMINES SHEET LENGTHS AND THICKNESSES USING THE COMPLETE 

* SOLUTION, INCLUDING GRAVITY, AND COMPARES THEM TO THE "nO 

* GRAVITY SOLUTION 


★ 




★ 

★ 

VARIABLE 

KEY 

* 

* 

WE 

: WEBER NUMBER 

•x 

* 

* 

FR 

: FRCUDE NUMBER 

* 

★ 

TAUB 

: TAU-BAR 

X 

* 

TIME3 

: TIME-3AR 

X 

★ 

Z3 

: Z-BAR 

X 

★ 

ZBNOG 

Z-3AR WITH NO GRAVITY 

f 

* 

NU 

: COEFFICIENT IN THE X-3AR EQUATION 

■x 


* ALL OTHER VARIABLES, TAUB ' S AND DIF'S, ARE JUST DUMMY VARIABLES * 

* IN THE ITERATIVE ROUTINE TO SOLVE FOR TAU-BAR * 

* DECLARE VARIABLES 

REAL WE, X3, DTAUB, TAU31 , TAUB 2 , DIF1, DIF2 , NU, TAUB, TAUBL, TAU3R, TAU3C, 
+DIFL, DIFR, DIFC, FR, TIMS3, Z3, Z3N0G 
INTEGER I,J 

* SET STEP SIZE FOR THE ROOT FINDER ROUTINE FOR TAUBAR * 

***'****************:r******************* i eir i eieie i r i r + if irie i r j [iril;i ' itiririfieiririfit j ritilrif ^' k 

DTAUB = 0.05 

*************************-k***l'**icic*******ieiri'**ir*irieir**it*ir**ir*-kir**iti'**-k-k** 

* SET HEADINGS FOR OUTPUT * 

a*****************************^*****************************^^^^^^^^^^^^ 

WRITE (8,*) ' DATA FOR SHEET LENGTHS AND THICKNESSES FOR WATER' 

WRITE { 8 , *) ' IN THE SMALL SCALE LSR. ' 

WRITE (8,*) ' ' 

WRITE (8,*) ' THE FOLLOWING PARAMETERS WERE USED:' 

WRITE ( 8 , * ) ' ' 

WRITE ( 8 , * ) ' SLIT WIDTH: 0.03416 M' 

WRITE(8,») ' SLIT THICKNESS: 8.128E-5 M' 

WRITE (8,*) ' DENSITY: 1000 KG/M^ ' 

WRITE ( 8 , * ) ' SURFACE TENSION: 7.282E-2 N/M' 

WRITE ( 8 , * ) • GRAVITY: 9.81 M/S~2' 

WRITE ( 8 , * ) ' ' 

WRITE (8,*) ' WEBER NUMBERS START AT 20 AND GO THROUGH 500* 

WRITE (8,*) ' ' 

WRITE (8,*) ' WEBER FROUDE X-BAR TAUBAR Z-BAR Z- 

+ BAR ' 

WRITE (8,*) ' NUMBER NUMBER NO 


m 


* . ’ T- „ * * 
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TAU31 =0.05 

30 TAU32 = TAUB1 + DTAU3 

DIF1 = 1.0 - NU* (TAU31** (-1.5) + 3 . 0 *TAUB1 * * { -0 . 5) - 4.0) - X3 
DIF2 - 1.0 - NU* (TAUB2** (-1.5) + 3 . 0 *TAUB2 * * (-0 . 5) - 4.0) - XB 
IF (DIF1 .EQ. 0.0) THEN 
TAUB = TAUB1 

ELSE IF (DIF2 .EQ. 0.0) THEN 
TAUB = TAUB 2 

ELSE IF (DIF1*DIF2 . LT . 0.0) THEN 
GOTO 40 

ELSE 

TAUB1 = TAUB2 
GOTO 30 
END IF 

40 TAU3L - TAUB1 
TAU8R = TAUB 2 

50 TAU3C - (TAUBL + TAUBR)/2.0 

DIFL - 1.0 - NU* (TAUBL** (-1.5) + 3 . 0* TAUBL* * (-0 .5) - 4.0) - X3 

DIFC - 1.0 - NU* (TAUBC** (-1.5) + 3 . 0 *TAUBC* * ( -0 . 5) - 4.0) - X3 

DIFR - 1.0 - NU* (TAUBR* * (-1.5) + 3 . 0*TAUBR* * (-0 .5) - 4.0) - XB 
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IF (ABS(DIFC) .LT. 0.0001) THEN 
TAUB - TAUBC 

ELSE IF (DIFL*D IFC . LT . 0.0) THEN 
TAUBR - TAUBC 
GOTO 50 

ELSE 

TAUBL = TAUBC 
GOTO 50 
END IF 

* SUBSTITUTE THE FOUND VALUE FOR TAU-BAR INTO THE TIME-BAR EQUATION * 

A*********************************************************************** 

TIMES = (FRMTAU3** (-1.0) - LOG (TAUB) - 1.0))/4.0 

★*********★*******★***★★****★**★*********★***★***★★★*******★**★****★**** 

* SUBSTITUTE THE VALUE FOR TIME-BAR INTO THE Z-BAR EQUATION * 

******»•». ***************************************************************** 

Z3 = (TIMEB*TIMEB) / (2 . 0*FR) + TIME3 

IF (X3 .EQ. 1.0) Z3 - 0.0 

************************************************************************ 

* OUTPUT THE DATA, END THE LOOPS, AND TERMINATE THE PROGRAM * 

*****»x»»*************************************************************** 

55 FORMAT (IX, F8.2,F10.3,F8.2,F11.5,F11.5,F11.5) 

WRITE (8, 55) WE, FR, X3, TAUB, ZB, ZBNCG 

20 CONTINUE 

WRITE {8,*) 1 * 

10 CONTINUE 

STOP 

END 
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***** + *** + '*'*************'* + ■+ + + + + + + *** + *'** + *** + + * + * + **'**************** + *** 

* AUTHOR: MATTHEW S. MCMASTSR * 

* DATS: 12-5-91 * 

* * 

* TITLE: SHAPE * 

* * 

* PURPOSE: TO DETERMINE THE CRCSS-SECTICNAL SHAPE CE THE SHEET * 

* ASSUMING PROPORTIONAL SECOND DERIVATIVES * 

* USING A FOURTH-ORDER RUNGS KUTTA TECHNIQUE * 

* * 

* VARIABLE KEY * 

* * 

* ETA : ETA ARRAY . * 

* THETA : DETA/DXI ARRAY ” * 

* XI : XI ARRAY * 

* AC3AR : AC -BAR * 

* Q : TAU/RE * 

* ALPHA : ALPHA * 

* ABAR : A- BAR * 

* ETAMAX : ETA-MAX * 

* ELLI? : ELLIPTICITY, 2 * ETA-MAX/ R-3AR * 

************************************************************************ 

************************************************************************ 

* DECLARE VARIABLES * 

ir*ir*'*-k m klririr'*-K-Kii:-kiric-kir*irie*1t1tit-K-*it-*iricit'xirir'k*'k-kic'*iririt-kirir*irie'*-x*'kic'k-k1c*1tic , k'k'k1t1c1rlt-x-x* 


' . ** 7 ' 

■ m 

\ i 

j m, 

•*; ; fe 

.*• % 
% 

i 


REAL ETA, THETA, XI, AC BAR, Q, ALPHA, ABAR, ETAMAX, DXI , DIF, FETA, FTHSTA, 
+F, K1 , X2 , K3, K4, Cl, C2 , C3 , C4, ELLI? 

INTEGER I,J 

DIMENSION ETA(0 : 6000) , THETA (0 : 6000 ) ,XI(0:6000) 

A*********************************************************************** 

* SET ACBAR., Q, ALPHA, AND THE XI ARRAY * 

**********-*»ir*************************'-*-************-** , **************'***-jr* 

ACBAR = 2000.0 

Q = 1 . 1 952 /ACBAR* *0 .3813 

ALPHA = -0.7571 + 0 . 120 6*LOG (ACBAR) 

DXI = 0.05 
DO 2, 1=0, 6000 
XI (I) = I*DXI 
2 CONTINUE 


************************************************************************ 

* USING THE RUNGE-KUTTA TECHNIQUE DETERMINE THE ETA AND THETA ARRAYS * 
************************************************************************ 


1 ETA(O) =0.0 

ETA ( I) - (DXI* ( (4.0* (1+ALPHA) ) /Q - DXI) ) **0.5 
THETA (0) =0.0 

THETA ( 1 ) = ( (2 . 0* (1 . 0+ ALPHA) ) /Q - DXI)/ETA(1) 

ABAR = (ETA (1) *DXI) /2 . 0 

DO 10, 1=2, 6000 
K1 = THETA(I-l) 

Cl = F (ETA ( I — 1 ) , THETA ( I — 1 ) , ABAR, ACBAR, ALPHA, Q) 
K2 = THETA { I — 1 ) + (DXI*Cl)/2.0 
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02 - F(ETA(I-1) + (DXI*Kl)/2.0, THETA (t-1) + (DXI *C1 ) /2 . 0 , ABAR, 
+ ACBAR, ALPHA, Q) 

K3 - THETA ( I — 1 ) + (DXI*C2)/2.0 

C3 - F (ETA(I-l) + (DXI*K2) /2 . 0, THETA ( 1-1 ) + (DXI *C2 ) /2 . 0 , AS AH, 

+ ACS AR, ALPHA, Q) 

K4 = THETA ( I-i) + (DXI*C3) 

C 4 = F (ETA (I-I) + (DXI *K3 ) , THETA ( 1-1 ) + (DXI *C3 ), ABAR, AC3AR, 

+ ALPHA, Q) 

ETA ( I ) =■ ETA(I-l) + DXI* (K1 + 2.0*K2 + 2.0*K3 + K4)/6.0 

THETA (I) =■ THETA (I-I) + DXI* (Cl + 2.0*C2 + 2.0*C3 + C4)/6.0 

IF ( ETA ( I ) .LT. ETA(I-l) .AND. ETA(I-2) .LT. ETA(I-l) .CR. ETA 
+ (I) .EQ. ETA(I-l) .AND. ETA.CI-2) . LT . ETA(I-l)) ETAMAX = ETA 

+ ( 1 - 1 ) 

IF (ETA (I) .GT. ETA ( I — 1 ) .AND. ETA(I-2) .GT. ETA(I-l) ,CR. ETA 
+ (I) .EQ. ETA(I-l) .AND. ETA(I-2) .GT. ETA(I-l)) GOTO 30 

ABAR = A3AR + DXI*(ETA(I) + ETA ( 1-1) ) /2 . 0 



10 CONTINUE 

* CHECK THE BOUNDARY CONDITION THAT AT THE MINIMUM ETA SHOULD EQUAL * 

* ONE AND CHANGE ALPHA UNTIL ETA DOES EQUAL ONE * 

* *★★***»********★*★******★*****★*★★*★★*****★★*■★**★*★*******★★★★*★■**★★★★*■ 

30 DIF = 1.0 - ETA(I-l) 

IF (ABS (DIF) .LT. 0.0001) GOTO 35 
ALPHA = ALPHA + DIF/20.0 
GOTO 1 

************************************************************************ 

* CHECK THE CONDITION THAT AT THE MINIMUM ABAR SHOULD EQUAL ACBAR AND * 

* CHANGE Q UNTIL ABAR DOES EQUAL ACBAR * 

************************************************************************ 

35 DIF = (ACBAR - ABAR) 

IF (ABS (DIF) /ACBAR .LT. 0.0001) GOTO 37 
Q =* Q + DIF/30000.0 
GOTO 1 

************************************************************************ 

* CALCULATE THE ELLIPTICITY, OUTPUT THE DATA, AND TERMINATE PROGRAM * 
************************************************************************ 

37 ELLIP =• (2.0*ETAMAX)/XI(I-1) 

38 FORMAT (1X,F8.2,F12.4,F12.4) 

WRITE (8,*) 'THE FOLLOWING ARE THE RESULTS FOR ETAMAX, RBAR, AND' 
WRITE (8,*) ' 2*ETAMAX/RBAR, GIVEN THE VALUE FOR ACBAR. ALSO' 

WRITE (8,*) 'GIVEN IS THE VALUES FOR TAU/R Z AND ALPHA NECESSARY’ 
WRITE (8,*) 'TO OBTAIN THE RESULT. ALSO GIVEN IS A TABLE OF' 

WRITE (8,*) 'SAMPLE VALUES OF ETA AND DETA/DXI FOR THE SHOWN' 

WRITE (8,*) 'VALUES OF XI.' 

WRITE (8,*) ’ ' 

40 FORMAT (1X,A,F8.3) 

41 FORMAT (IX, A, F9 . 4) 

42 FORMAT ( IX, A, F10 . 5) 

WRITE (8,40) 'ACBAR - ', ACBAR 



. *■ V v V I* 

v. - 


22 


WRITS 

(8,40) 

'ETAMAX ETAMAX 

WRITE 

(8,40) 

* RBAR - ' , XI (1-1) 

WRITE 

(8,41) 

' 2 ’ ETAMAX /RBAR - ' , ELLIP 

WRITE 

(8,42) 

'TAU/RE -',Q 

WRITE 

(8,42) 

■ALPHA ALPHA 

WRITE 

(8,*) » 

» 

WRITE 

(8,*) ' 

XI ETA 

WRITE 

(8,*) » 

— 

DO 22, 

J - 0, 

1-1,20 


WRITE (8,38) XI (J) , ETA (J) , THETA (J) 
CONTINUE 



DETA/DXT ’ 


STOP 

END 


* ♦★★★★★★♦★*****************ir**********************'*'' r *******************^ 

* FUNCTION F SUBPROGRAM * 

* F IS SASICLY THE EQUATION BEING SOLVED * 

*************************■*■*****»***■****************'*****************»»»» 


FUNCTION F (FETA, FTHETA, ABAR, AC3AR, ALPHA, Q) 

REAL FETA, FTHETA, ABAR, ALPHA, Q, F, AC3AR 

F = ( ( ( (ABAR*ABAR* (1.0 + FTHETA*FTHETA) ) / (FETA*FETA*AC3AR 
+ *AC3AR) - Q) *(1.0 + FTHETA*FTHSTA) **1.5) / (2 .0* (1.0 + ALPHA))) 

END 


X 
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a**************************************************************** 

* AUTHOR: MATTHEW S. MCMASTER 

* DATE: 12-6-91 




* TITLE: SHAPE 2 * 

* * 

* PURPOSE: TO DETERMINE THE CROSS-SECTIONAL SHAPE OF THE SHEET * 

* ASSUMING A CONSTANT SECOND DERIVATIVE * 

* USING A FOURTH ORDER RUNGE-KUTTA TECHNIQUE * 

* • 

* variable: key * 


* ETA 

* THETA 

* XI 

* ACBAR 

* Q 

* R 

* A3AA 

* ETAMAX 

* ELLIP 


ETA ARRAY 
DETA/DXI ARAY 
XI ARRAY 
AC -BAR 
TAU/RE 

D A 2ETA/DZETA A 2 

A-3AR 

ETA-MAX 

ELLIPTICITY, 2* ETA-MAX/ R-3AR 


★ ★★***★*»»*★★**★*★★****★*»*********★★***★**★★******★******★★**'★****'*'*'**'* 
* declare: variables * 

**ir**x*'x'r‘*-xx*x**it1tx**ic***-)rir + *ic******itie*ieie********'*-r* + **********ir*'**x**xx 

REA. L ETA, THETA, A3AR, Q, ACBAR, R, DXI, DIF, ETAMAX, F, XI, FETA, FTHSTA, 

+K1, K2,K3,K4, Cl, C2,C3,C4, ELLIP 
INTEGER I, J 

DIMENSION ETA(0 : 60CC) , THETA ( 0 : 6000 ) ,XI(0:6000) 


* SET ACBAR, Q, R, AND THE XI ARRAY , * 

*★*****★*»*********★★★****»*★★*■*★★'*■**★★*■★********★★***★★★*★*★**'★*★■«■■***** 


ACBAR = 2000.0 

Q = 403. 65 06/ ACBAR* *1. 3278 

R = -0.4827 + 0.0791*LCG(AC3AR) 

DXI =■ 0.05 
DO 2, I - 0,6000 
XI (I) - I *DXI 
2 CONTINUE 

************************************************************************ 
* USING THE RUNGE-KUTTA TECHNIQUE DETERMINE THE ETA AND THETA ARRAYS * 

********************************************************************** 1 ,* 

1 ETA(O) - 0.0 

ETA { 1) - (DXI* (4 . 0/Q - DXI) ) **0 . 5 
THETA (0) -0.0 

THETA ( 1 ) - (2. 0/Q - DXI) /ETA ( 1) 

ABAR - (ETA(l) *DXI)/2.0 

DO 10, I - 2,6000 
K1 - THETA ( 1-1) 

Cl - F (ETA ( 1-1) , THETA (1-1) , ABAR, ACBAR, R,Q) 

K2 - THETA ( I — 1 ) + (DXI*Cl)/2.0 


. . ..v * .i -."L. y ■ ** tXZ * " 




C2 

- F(ETA(I-1) 

+ 

(2X1 *K1) /2. 0, THETA (I- 


AC3AR, R, 

Q) 


K3 

- THETA ( X — 1 ) 

+ 

(2XI*C2) /2.0 

C3 

- F (ETA { 1-1) 

+ 

(2X1 *K2) /2 . 0, THETA ( I- 


ACBAR, R, 

Q) 


K4 

= THETA (1-1) 


(2X1 *C3) 

Cl 

- F (ETA ( 1-1) 


(2X1 *K3) , THETA (1-1) + 


R,Q) 

ETA(I) - ETA(I-l) + 2X1* (Kl + 2.0*K2 + 2.‘. 
THETA (I) - THETA(I-l) + DXI*(C1 + 2.0*C2 - 


+ 

+ 

+ 


IF (ETA (I) .LT. ETA(r-l) .AND. ETA (1-2) 
(I) .EQ. ETA ( I — 1 ) . AND . ETA(I-2) .LT. E 
( 1 - 1 ) 

IF (ETA (I) .GT. STA(I-l) .AND. ETA (1-2) 
(I) .EQ. ETA(I-i) .AND . ETA (1-2) .GT. E' 



A3AR - ABAR + DXI* (ETA (I) + ETA(I-l) ) / 2 . D 


10 CONTINUE 


* CHECK THE BOUNDARY CONDITION THAT AT THE MINIMUM EH FEU 

* CNE AND CHANGE R UNTIL ETA. ICES EQUA.L CNS 

ic'k-k-k'X'r-klt-kir'k-kltitltiC'Xi'icic'k-x-xicic'icicTr-K'Kv-K-K'Kici'irir'x-K-xw-x'Kitiriririt-K'K'w: ti t » tt t ■ 


30 DIF = 1.0 - ETA(I-l) 

IF (ABS(DIF) .LT. 0.0CCI) GOTO 35 
R = R - DIF/30.0 
GOTO 1 

icitit'k'K-X'K**ic-kir'k'kiC'kit'*'k-kitic'Xifk* + -k-wK-x-K-x'Xit'kit'*it'K-**Tir*itit*ir** tT-i »m t 

* CHECK THE CONDITION THAT AT THE MINIMUM A3AR SHOULD EQUAL 

* CHANGE Q UNTIL ABAR DOES EQUAL AC3AR 

+ -»****■*■'*'«■■*▼ t-tt t-t ▼ t 

35 DIF = AC3AR - ABAR. 

IF (A3S (DIF) /AC3.AR .LT. 0.0001) GOTO 37 
Q = Q - DIF/35000.0 
GOTO 1 


* CALCULATE THE ELLIPTICITY, CCTPUT THE DATA, AND TEP2 
*★★**•**■■*★********★**■**■**★*★*■*»*»-**********■-*•*»★***■** 


37 ELLIP - (2 . 0*ETAMAX) /XI (1-1) 

38 FORMAT (IX, F8 . 2, F12 . 4, F12. 4) 

WRITS (8,*) 'THE FOLLOWING ARE THE RESULTS FOR STSfP.T, 
WRITE (8,*) '2*ETAMAX/S=AE, GIVEN THE VALUE FOR ADE3A- 
WRITE (8,*) 'GIVEN IS THE VALUES FOR TAU/RE AND ITIETA 
WRITE (8,*) 'NECESSARY TO OBTAIN THE RESULT. ALSO HITE 
WRITE (8,*) 'TABLE OF SAMPLE VALUES OF ETA AND DS3.UXI 
WRITE (8,*) 'SHOWN VALUES OF XI.' 

WRITE (8,*) * ' 

40 FORMAT (1X,A,F8.3) 

41 FORMAT (IX, A, F9 . 4) 

42 FCPMAT (IX, A,F10.5) 

WRITE (8,40) ' ACBAR - \AC3AR 
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WRITE (8,40) 'ETAMAX - ' , ETAMAX 
WRITS (8,40) 'REAR -',XI(I-1) 

WRITE (8,41) *2 ‘ETAMAX/ RBAR - ’ , ELLIP 

WRITE (8,42) 'TAU/RE - ' , Q 
WRITS (8,42) 'D~2ETA/DZETA"2 ** ' , R 
WRITE (8,*) ' » 

WRITE (8,*) ' XI ETA DETA/DXI' 

WRITS (8,*) • — , 

DO 22, J - 0,1-1,20 

WRITE (8,38) XI (J) , ETA (J) , THETA (J) 

22 CONTINUE 

STOP 

END 


************************************************************************ 

* FUNCTION F SUBPROGRAM * 

* F IS 3ASICLY THE EQUATION BEING SOLVED * 

*********** *************************************************»******, ^.^^^ 

FUNCTION F (FETA, FTHETA, ABAR, AC3AR, R, Q) 

REAL FETA, FTHETA, ABAR, R,Q,F,ACBAR 

- = ( ( (ABAR*ABAR‘ (1 .0 + FTHSTA*FTHSTA) ) / (FETA*FETA*AC3AR* 

+ AC3AR) - Q) *(1.0 + FTHETA*FTHETA) **1.5) /2 . 0 - R 


END 


M 








ct 
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